Outubro de 2022
Uma série temporal (ST) é uma coleção de observações feitas sequencialmente no tempo.
Muitas STs que conhecemos estão no campo da economia, as cotações diárias do mercado de ações, índice de inflação ou taxa de desemprego mensal.
Em ciências sociais seguem séries populacionais, como taxas de natalidade.
Na computação, a evolução de um software pode ser estudado através de suas releases.
Na epidemiologia surge o número de casos de dengue ao longo de um período de tempo.
Na medicina, a ressonância magnética de ondas cerebrais pode ser usada para estudar como o cérebro reage a certos estímulos experimentais.
Na área de saúde, o comportamento de uma série temporal ficou mais conhecido com os casos de COVID-19 apresentados diariamente.
Figura 1: Número mensal de manchas solares.
R, os comandos são:sunspots.m <- data.frame(year = index(sunspot.month),
value = melt(sunspot.month)$value)
ggplot(sunspots.m, aes(x = year, y = value)) +
geom_line(colour = "red")+
ylab("Número de manchas solares") +
xlab("Ano")
A principal diferença entre séries temporais e dados transversais está na estrutura de dependência.
A análise de dados transversais trata, principalmente, observações i.i.d – independentes e identicamente distribuídas.
Enquanto que em séries temporais, cada nova observação é estocasticamente dependente de observações anteriores.
Ao mesmo tempo que a dependência é nossa maior inimiga, ela será nossa grande aliada.
Em estudos de séries temporais é possível avaliar, através de medições de temperaturas, a hipótese de aquecimento global.
É possível avaliar também se os níveis de poluição podem influenciar a mortalidade diária em grandes metrópolis.
Os dois últimos são exemplos de dados experimentais que podem ser usados para ilustrar o processo pelo qual a metodologia estatística clássica é aplicada em séries temporais correlacionadas.
De modo geral, o primeiro passo em qualquer investigação de séries temporais é a observação dos dados traçados ao longo do tempo.
É importante mencionar que existem duas abordagens distintas, mas não mutuamente excludentes, na análise de séries temporais:
A abordagem no domínio do tempo é motivada pela presunção de que a correlação entre pontos adjacentes no tempo é melhor explicada em termos de uma dependência entre o valor atual e valores passados.
Assim, busca-se modelar valores futuros de uma série temporal como uma função paramétrica dos seus valores atuais e passados.
Neste cenário, podemos usar regressões lineares do valor presente de uma série temporal sobre seus próprios valores passados e sobre valores passados de outras séries.
Por outro lado, a abordagem do domínio da frequência pressupõe que as características de interesse em análises de séries temporais se relacionam com variações sistemáticas encontradas naturalmente na maioria dos dados.
Essas variações periódicas são frequentemente causadas por fenômenos biológicos, físicos ou ambientais:
Os estudos da periodicidade são comuns na economia e nas ciências sociais, onde há interesse na periodicidade anuais de séries como o desemprego mensal ou na taxa de natalidade mensal.
Alguns dos problemas e questões de interesse do analista de séries temporais podem ser apresentados através de dados reais extraídos de diferentes áreas de concentração.
Os exemplos a seguir ilustram alguns tipos de dados de séries temporais, bem como algumas questões estatísticas que podem surgir sobre esses dados.
Esta série temporal refere-se aos dados de vendas mensais de medicamentos antidiabéticos na Austrália de 1992 a 2008.
São 204 meses. Total de scripts mensais para produtos farmacêuticos abrangidos pelo código A10 do ATC, tal como registado pela Australian Health Insurance Commission.
A série a10 pertence ao pacote fpp.
Figura 2: Vendas mensais de medicamentos antidiabéticos.
R, os comandos são:ggplot(a10, aes(as.Date(a10), as.matrix(a10))) +
geom_line(colour = "red")+
ylab("Vendas") +
xlab("Tempo")
Esta série temporal refere-se aos dados mensais de mortes por pneumonia e gripe a cada 10.000 pessoas nos Estados Unidos por 11 anos, de 1968 a 1978.
São 132 meses.
A série flu pertence ao pacote astsa.
Figura 3: Mortes por pneumonia e gripe a cada 10.000 pessoas.
ggplot(flu, aes(as.Date(flu), as.matrix(flu))) +
geom_line(colour = "red")+
ylab("Taxa") +
xlab("Tempo")
Esta série temporal refere-se aos dados casos de AIDS notificados no SINAN, declarados no SIM e registrados no SISCEL/SICLOM(1), segundo capital de residência por ano de diagnóstico. Brasil, 1980-2021.
Fonte: MS/SVS/Departamento de Doenças de Condições Crônicas e Infecções Sexualmente Transmissíveis (DCCI).
> AIDS <- read.csv("Dados SUS/AIDS.csv", sep = ";")
> AIDS <- ts(AIDS$Frequência, start=1980, frequency=1)
São 41 anos.
A série foi coletada no site do DATASUS.
Figura 4: Casos de AIDS notificados no SINAN.
R, os comandos são:ggplot(AIDS, aes(as.Date(AIDS), as.matrix(AIDS))) +
geom_line(colour = "red")+
ylab("Casos") +
xlab("Tempo") +
scale_y_continuous(labels = addUnits)
scale_y_continuous(labels = addUnits) utilizada no gráfico anterior roda com base na função abaixoaddUnits <- function(n) {
labels <- ifelse(n<1000,n,#less than thousands
ifelse(n<1e6,paste0(round(n/1e3),'k'),#in thousands
ifelse(n<1e9,paste0(round(n/1e6),'M'),#in millions
ifelse(n<1e12,paste0(round(n/1e9),'B'),#in billions
ifelse(n<1e15,paste0(round(n/1e12),'T'),#in trillions
'too big!'
)))))
return(labels)
}
Podemos ter interesse em analisar várias séries temporais ao mesmo tempo.
Considere as séries de um estudo para investigar possíveis associações entre mortalidade e poluição no condado de Los Angeles.
Neste exemplo, vamos considerar as séries temporais no período de 1970 a 1979:
Ambas as séries com 508 semanas. Para maiores detalhes http://www.sungpark.net/ShumwayAzariPawitan88.pdf
cmort e part pertencem ao pacote astsa.Figura 5: Mortalidade cardiovascular.
autoplot(cbind(Mortality=cmort), xlab='Tempo',
ylab='Mortalidade', col=2)
autoplot(cbind(Mortalidade=cmort, Partículas=part),
xlab='Tempo', ylab='')
Ambas as séries tendem a apresentar comportamento repetitivo, com ciclos regulares.
Parece haver uma relação entre as duas séries.
Esta possibilidade sugere o uso de técnicas da análise de regressão como um procedimento para relacionar as duas séries.
Dado que temos uma ou mais séries temporais para realizar a análise, como podemos fazê-la?
A característica peculiar de dados de séries temporais é que as observações não são independentes.
Logo, qualquer análise deve considerar a ordem em que as observações são coletadas.
Cada observação numa série é uma observação bivariada com o tempo sendo a segunda variável.
Descrição. Descrever os dados usando estatísticas descritivas e métodos gráficos;
Modelagem. Encontrar um modelo estatístico adequado para descrever o processo gerador de dados. Obviamente, todos os modelos são aproximações e a construção de modelos é uma arte;
Previsão. Estimar (prever) valores futuros da série.
Em análise de séries temporais é crucial desenvolver modelos estatísticos que forneçam descrições plausíveis para os dados.
Encontrar um modelo estatístico capaz de descrever características dos dados que aparentemente flutuam de forma aleatória ao longo do tempo.
Uma série temporal pode ser definida como uma coleção de variáveis de acordo com a ordem em que são obtidas no tempo.
Por exemplo, uma série temporal é uma sequência de variáveis, \(x_1, x_2,\dots\), em que \(x_1\) denota o valor observado no primeiro ponto de tempo, assim por diante.
Consideremos aqui a série temporal \(\{x_t\}\), indexada por \(t = 1, 2,\dots\), tipicamente discreto.
Podemos imaginar que, em cada uma das séries apresentadas antes, há a suposição de que pontos adjacentes no tempo estão correlacionados.
Ou seja, o valor no tempo \(t\), digamos \(x_t\), depende dos valores passados \(x_{t-1}, x_{t-2},\dots\)
Para ilustrar, vamos considerar o exemplo a seguir.
\[w_t \sim wn(0, \sigma_w^2).\]
\[w_t \sim {\rm N}(0, \sigma_w^2).\]
set.seed(1) wn <- rnorm(500,0,1) head(wn)
## [1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078 -0.8204684
ts_wn <- data.frame(Time = seq(0, length(wn)-1, 1), WN = wn) head(ts_wn) #mostra as 6 primeiras observações
## Time WN ## 1 0 -0.6264538 ## 2 1 0.1836433 ## 3 2 -0.8356286 ## 4 3 1.5952808 ## 5 4 0.3295078 ## 6 5 -0.8204684
tail(ts_wn)#mostra as 6 últimas observações
## Time WN ## 495 494 3.8102767 ## 496 495 -1.1089100 ## 497 496 0.3075666 ## 498 497 -1.1068945 ## 499 498 0.3476536 ## 500 499 -0.8732645
Figura 6: Série ruído branco gaussiano.
Lembram qual o pressuposto sobre os resíduos nos modelos clássicos de regressão?
Portanto, devemos considerar a correlação serial.
Duas formas de introduzir a correlação serial em modelos de séries temporais são exemplificadas a seguir.
Podemos suavisar a série de ruído branco \(w_t\) através do modelo média móvel.
Por exemplo, considere substituir \(w_t\) no exemplo anterior por uma média de seu valor atual e seus vizinhos imediatos no passado e no futuro.
\[v_t = \frac{w_{t-1} + w_t + w_{t+1}}{3}\]
Figura 7: Média móvel de três pontos da série ruído branco gaussiano.
gg_wn <- ggplot(data = ts_wn, aes(x = Time, y = WN)) +
geom_line() + labs(x = "\n Tempo ", y = "RB \n") +
theme(plot.title = element_text(face="bold",
colour="darkblue",
size = 12,hjust = 0.5))
gg_wn
ts_mm <- data.frame(Time = seq(0, length(wn)-1, 1),
MM = rollapply(wn,3,mean,align='right',
fill=NA))
gg_mm <- ggplot(data = ts_mm, aes(x = Time, y = MM)) +
geom_line() + labs(x = "\n Tempo ", y = "MM \n") +
theme(plot.title = element_text(face="bold",
colour="darkblue",
size = 12,hjust = 0.5))
gg_mm
grid.arrange(gg_wn,gg_mm, ncol = 1)
Examinando o gráfico, vemos que \(v_t\) é uma versão mais suavizada com relação à primeira série \(w_t\).
Ou seja, na média móvel as oscilações são mais lentas e algumas das oscilações mais rápidas são extraídas.
\[x_t = x_{t-1} - 0.9 x_{t-2} + w_t,\] para \(t=1, 2, \dots, 500\).
A equação acima representa a regressão do valor atual \(x_t\) de uma série temporal em função dos dois últimos valores da própria série; por isso, o termo autorregressão.
A série resultante desta autorregressão é mostrada na figura a seguir.
Nota-se alguma periodicidade na série.
O modelo autorregressivo é bastante utilizado no ajuste de muitas séries observadas e serão estudados em detalhes mais adiante.
Para gerar (simular) e plotar estes dados autorregressivos no R:
set.seed(2012)
yt <- arima.sim(list(order=c(2,0,0), ar=c(1,-.9)), n=500)
ts_ar <- data.frame(Time = seq(0, length(yt)-1, 1), yt)
gg_ar <- ggplot(data = ts_ar, aes(x = Time, y = yt)) +
geom_line() +
labs(x = "\n Tempo ", y = "AR \n", title = "Autoregressão \n") +
theme(plot.title = element_text(face="bold", colour="darkblue",
size = 12,hjust = 0.5))
gg_ar
\[ \rho(s,t) = \frac{\gamma(s,t)}{\sqrt{\gamma(s,s) \gamma(t,t)}}. \]
O ACF mede a previsibilidade linear da série no tempo \(t\), digamos \(x_t\), usando apenas o valor \(x_s\).
Tem-se que \(-1 \leqslant \rho(s,t) \leqslant 1\).
Se é possível prever perfeitamente \(x_t\) a partir de \(x_s\) através da relação linear \(x_t = \beta_0 + \beta_1 x_s\), então:
\(\rho(s,t) = +1\), quanto \(\beta_1 > 0\);
\(\rho(s,t) = -1\), quanto \(\beta_1 < 0\);
Então, \(\rho(s,t)\) é uma medida aproximada da capacidade de prever a série no instante \(t\) com o valor no instante \(s\).
\[ \widehat\rho(h) = \frac{\sum\limits_{t=1}^{n-h} (x_{t+h} - \bar{x})(x_t - \bar{x})}{\sum\limits_{t=1}^{n} (x_t - \bar{x})^2}, \] para \(h = 0, 1,\dots, n-1\).
Figura 8: Dispersão da série cmort e defasagens; \(h=1\) e \(h=6\).
A Figura 8 mostra um exemplo usando a série de Mortalidade Cardiovascular (cmort) onde \(\widehat\rho(1) = 0.772\) e \(\widehat\rho(6) = 0.528\).
O código a seguir foi usado para construir a Figura 8.
rho_ = round(acf(cmort, 6, plot=FALSE)$acf[-1], 3) # valores da ACF cmort.lag1 <- dplyr::lag(as.matrix(cmort),1) cmort.lag6 <- dplyr::lag(as.matrix(cmort),6) df_cor <- data.frame(cmort.lag1, cmort.lag6, as.matrix(cmort)) g1<-ggplot(data = df_cor, aes(cmort.lag1, cmort)) + geom_point(color='blue') + geom_label(data = NULL, x = 80, y = 125, label = rho_[1], size=5, col="red") g2<-ggplot(data = df_cor, aes(cmort.lag6, cmort)) + geom_point(color='blue') + geom_label(data = NULL, x = 80, y = 125, label = rho_[6], size=5, col="red") grid.arrange(g1,g2, ncol=2, nrow=1)
Propriedade: Se \(x_t\) é ruído branco, para \(n\) grande, a ACF amostral, \(\widehat\rho_x(h)\), para \(h = 1, 2,\dots, H\), é aproximadamente normal com média zero e desvio-padrão dado por \(\frac{1}{\sqrt{n}}\).
Importante: É devido a esta propriedade que se aplicam técnicas para transfomar uma série temporal em uma série de ruído branco, de modo a não violar a suposição de alguns modelos estatísticos.
Figura 9: ACF da série Mortalidade Cardiovascular (cmort)
Figura 10: ACF da série Partículas de Poluição (part).
A ACF de cmort exibe periodicidades correspondentes a valores separados por 17 unidades de tempo, que correspondem a 17 semanas ou 1 quadrimestre.
Observações das semanas referentes aos \(1^{o}\) e \(3^o\)s quadrimestres de cada ano estão correlacionadas positivamente;
Observações das semanas referentes aos \(2^o\)s de cada ano estão negativamente correlacionadas.
A ACF de part exibe periodicidades correspondentes a valores separados por 12 unidades de tempo, que correspondem a 12 semanas ou 1 trimestre.
ggAcf(cmort, 104, main="Mortalidade Cardiovascular") ggAcf(part, 104, main="Partículas de Poluição")
Até aqui, não fizemos qualquer suposição especial sobre o comportamento da série temporal.
Os exemplos anteriores sugerem alguma regularidade ao longo do tempo no comportamento da série.
Assim como nos modelos clássicos de regressão, quando pretendemos utilizar modelos para descrever séries temporais, é necessário introduzir suposições simplificadoras.
A principal suposição é a de estacionariedade.
Intuitivamente, uma série \(x_t\) é estacionária se ela se desenvolve no tempo de modo que a escolha de um instante no tempo não é importante (Morettin e Toloi, 2006).
Em outras palavras, as características de \(x_{t+s}\), para todo \(s\), são as mesmas de \(x_t\).
Tecnicamente, há duas formas de estacionariedade:
Para mais detalhes, sugiro ver Morettin e Toloi (2006).
O modelo de regressão linear é mais utilizado no contexto de dados clássicos do que em séries temporais.
Estes modelos de regressão são importantes para compreender os modelos de séries temporais que discutiremos mais adiante.
Já vimos que os modelos autorregressivos expressam \(x_t\) como uma combinação linear de valores passados \(x_{t-1}, x_{t-2},\dots, x_{t-p}\) da própria série.
Além disso, \(x_t\) pode depender de valores defasados de outra série, digamos \(y_{t-1}, y_{t-2},\dots, y_{t-q}\), que a influencia.
Assumimos que a série temporal dependente \(x_t\), para \(t = 1,\dots, n\), é influenciada por um conjunto de séries independentes, \(z_{t1}, z_{t2},\dots, z_{tq}\), conhecidas.
Esta é uma condição necessária para aplicar a regressão linear convencional.
Então, expressamos essa relação através do modelo de regressão linear
\[ x_t = \beta_1 z_{t1} + \beta_2 z_{t2} + \cdots + \beta_q z_{tq} + w_t, \] em que \(\beta_1, \beta_2,\dots, \beta_q\) são os coeficientes desconhecidos a serem estimados, \(w_t\) é o erro aleatório (ou ruído) independente e normalmente distribuído com média 0 e variância \(\sigma_w^2\).
Lembrem que nesta série há uma aparente tendência ascendente que é usada para argumentar a hipótese do aumento dos casos de AIDS.
Podemos usar regressão linear para estimar essa tendência ajustando o modelo simples
\[ x_t = \beta_1 + \beta_2 t + w_t,\qquad t = 1980, \dots, 2020. \]
Este é um modelo de regressão, onde \(t\) é a variável independente.
Observe que estamos fazendo a suposição de que os erros, \(w_t\), são uma sequência iid normal.
# regressão dos casos de AIDS no tempo modelo = lm(AIDS~time(AIDS)) summary(modelo)$coefficients
## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -2205365.209 233235.3452 -9.455536 1.218968e-11 ## time(AIDS) 1115.431 116.6156 9.565020 8.878637e-12
\[\widehat{x}_t = -2205365.21 + 1115.431\ t.\]
Figura 11: Casos de AIDS com tendência linear estimada.
Observem que a linha de tendência estimada via regressão linear não capta bem a tendência dos dados.
Abaixo temos os comandos usados para contruir a Figura 11.
ggplot(data = AIDS, aes(as.Date(AIDS), as.matrix(AIDS))) + geom_point(color='blue') + labs(x = "\n Tempo ", y = "Casos \n") + scale_y_continuous(labels = addUnits)+ geom_smooth(method = "lm", se = FALSE, color='red')
Em geral, é desejável que a série temporal seja estacionária, de modo que a sua média não mude com o tempo.
Além disso, a dependência entre os valores da série é algo que devemos medir.
No entanto, é difícil medir esta dependência se a estrutura dela estiver mudando com o tempo.
Então, para obter qualquer análise estatística significativa da série, é crucial que a média e a função de autocovariância satisfaçam as condições de estacionaridade:
a média, \(\mu_t\), é constante e não depende do tempo \(t\), e
a função de autocovariância, \(\gamma(s,t)\), depende de \(s\) e \(t\) somente através de sua diferença \(\vert s-t \vert\).
Muitas vezes, este não é o caso, e vamos mencionar alguns métodos nesta seção para minimizar os efeitos da não-estacionaridade.
Uma forma comum de lidar com a não-estacionaridade é quando o processo tem comportamento estacionário em torno de uma tendência.
Podemos escrever este tipo de modelo como
\[x_t = \mu_t + y_t,\] em que \(x_t\) é a série, \(\mu_t\) denota a tendência e \(y_t\) é um processo estacionário.
Logo, pode ser vantajoso remover a tendência como um primeiro passo em uma análise exploratória de séries temporais.
Para isso, basta obter uma estimativa razoável do componente de tendência, digamos \(\widehat{\mu}_t\), e então trabalhar com os resíduos
\[ \widehat{y}_t = x_t -\widehat{\mu}_t. \]
\[ \widehat{\mu}_t = -2205365.21 + 1115.431\ t. \]
\[ \begin{aligned} \widehat{y}_t &= x_t -\widehat{\mu}_t\\[1mm] \widehat{y}_t &= x_t + 2205365.21 -1115.431\ t. \end{aligned} \]
Figura 12: Eliminando tendência da série dos casos de AIDS.
modelo = lm(AIDS~time(AIDS), na.action=NULL) residuo = resid(modelo) ggplot(data = residuo, aes(as.Date(residuo), as.matrix(residuo))) + geom_line(color='blue') + labs(x = "\n Data ", y = "Resíduos \n")
R, plotar a série sem tendência e as ACF amostrais anteriores, usamos os seguintes comandos:modelo = lm(AIDS ~ time(AIDS), na.action=NULL) ggplot(data = AIDS, aes(as.Date(AIDS), as.matrix(AIDS))) + geom_point(color='blue') + labs(x = "\n Data ", y = "Casos \n") plot(resid(modelo), type="o", main="sem tendência") par(mfrow=c(3,1)) # plot ACFs ggAcf(AIDS, 48, main="Casos de AIDS") ggAcf(resid(modelo), 48, main="sem tendência")
Outra prática comum para tornar uma série estacionária é tomar as diferenças sucessivas.
Uma vantagem de diferenciar a série sobre o processo de remover tendência anterior é que nenhum parâmetro é estimado na operação de diferenciação.
A primeira diferença é denotada por
\[ \nabla x_t = x_t - x_{t-1}. \]
\[ \begin{aligned} \nabla^2 x_t &= \nabla(\nabla x_t)\\[1mm] &= \nabla(x_t - x_{t-1})\\[1mm] &= (x_t - x_{t-1}) - (x_{t-1} - x_{t-2})\\[1mm] &= x_t - 2x_{t-1} + x_{t-2}, \end{aligned} \] pode eliminar uma tendência quadrática, e assim por diante.
\[ \nabla^d x_t = \nabla(\nabla^{d-1} x_t). \]
Em situações normais, será suficiente tomar uma ou duas diferenças para que a série se torne estacionária.
A técnica de diferenciação é um componente importante do modelo ARIMA discutido mais adiante.
autoplot(diff(AIDS), type="o", main="primeira diferença")
Parece que a série da primeira diferença já mostra autocorrelação mínima.
Para obter e plotar os dados da primeira diferença, e gerar a ACF amostral, usamos os comandos:
plot(diff(AIDS), type="o", main="primeira diferença") ggAcf(diff(AIDS), 48, main="primeira diferença")
Esta é outra técnica exploratória de dados usada para visualizar relações entre diferentes defasagens de séries temporais.
Anteriormente, vimos que a ACF nos indica existência de relação linear entre a série e seus próprios valores defasados.
No entanto, a ACF pode mascarar uma possível relação não-linear entre \(x_t\) e seu valores passados, \(x_{t-h}\).
Para avaliar possíveis relações não-lineares via matriz de dispersão, podemos plotar a série original, \(x_t\), contra suas defasagens \(x_{t-h}\).
# install.packages("astsa")
library(astsa)
lag1.plot(cmort, 12)
Figura 13: Matriz de dispersão de valores atuais de cmort, \(M_t\), contra defasagens \(M_{t-h}\).
As autocorrelações amostrais são exibidas no canto superior direito nos gráficos de dispersão.
Os gráficos também exibem as linhas de alisamento de dispersão (lowess) usadas para indicar a relação não-linear.
Observação: Discutiremos alisamento a seguir, mas entenda lowess como um método para ajustar a regressão local.
Além disso, vemos fortes relações lineares positivas entre \(M_t\) e as defasagens \(M_{t-1}, M_{t-2}, M_{t-3}, M_{t-4}\) e depois começa a reduzir o valor da correlação à medida que aumenta a quantidade de defasagens.
Da mesma forma, podemos avaliar possíveis relações não-lineares entre a série de mortalidade cardiovascular, \(M_t\), e desfasagens do nº de partículas, \(P_{t-h}\).
Lembrem do exemplo onde queríamos prever a média de mortalidade semanal a partir da quantidade de partículas de poluição.
Figura 14: Matriz de dispersão de valores atuais da mortalidade cardiovascular, \(M_t\), contra defasagens de partículas, \(P_{t-h}\).
A figura anterior mostra a matrix de dispersão da série de mortalidade cardiovascular contra a série do nº de partículas e suas defasagens, \(P_{t-h}\).
Além disso, a figura exibe as correlações cruzadas, bem como os ajustes lowess.
Nota-se relação não-linear de forma moderada entre mortalidade cardiovascular \(M_t\) e partículas de poluição nas defasasgens \(P_{t-4}, P_{t-6}, P_{t-7}, P_{t-8}\), indicando que a série de partículas pode influenciar antecipadamente a série de mortalidade cardiovascular.
astsa e o código no R éinstall.packages("astsa")
library(astsa)
astsa::lag1.plot(part, 12) # dispersão com uma série
astsa::lag2.plot(part, cmort, 8) # dispersão com duas séries
Este método é útil para descobrir algumas características em séries temporais, como tendências de longo prazo e componentes sazonais.
Em particular, se \(x_t\) representa as observações, então \[ m_t = \sum\limits_{j=-k}^{k} a_j x_{t-j}, \] em que \(a_j = a_{-j} \geqslant 0\) e \(\sum_{j=-k}^{k} a_j = 1\) é a média móvel dos dados.
Figura 15: Óbitos diários por COVID-19 na Bahia
#library(readxl)
covid_bahia <-
read_excel("Dados/Banco Estadual COVID-19 ÓBITOS_17-10-2022.xlsx")
covid_bahia$Mortes <- seq(1,1,length(covid_bahia$`DATA OBITO`))
#View(covid_bahia)
class(covid_bahia$`DATA OBITO`)
### Ajustando a variável data usando
## o pacote lubridate
## Indicar a ordem da data em inglês
## Y-year , m- month d-day
#library(lubridate)
covid_bahia$`DATA OBITO` <-
as.Date(parse_date_time(covid_bahia$`DATA OBITO`,
"%y/%m/%d"))
class(covid_bahia$`DATA OBITO`)
df_covid = tibble(covid_bahia)
dplyr::glimpse(df_covid)
df_covid <- df_covid%>%
dplyr::group_by(`DATA OBITO`)%>%
dplyr::summarise_if(is.numeric,sum)
#library(zoo)
df_covid=df_covid %>%
mutate('Média Móvel'=rollapply(Mortes,14,
mean,align='right',
fill=NA))
## Agrupando para termos categorias de
## Médias Móvel e Óbitos no mesmo gráfico
### Deixando em outro conjunto Bahia
df_covid2=df_covid %>%
gather(c('Média Móvel','Mortes'),key="Séries",
value="Valor")
ggplot(df_covid2,aes(x=`DATA OBITO`,y=Valor,fill=Séries,
colour=Séries))+
geom_line()+
labs(x="",y="Número de óbitos")
Outra abordagem para suavizar uma série é a regressão do vizinho mais próximo.
A técnica baseia-se na regressão de \(k\)-vizinhos mais próximos, onde se usa apenas os dados \(\{x_{t-k/2},\dots, x_t,\dots, x_{t+k/2}\}\) para prever \(x_t\) via regressão, e então define \(m_t = \widehat{x}_t\).
Lowess é um método de alisamento complexo, mas a ideia básica é próxima à regressão do vizinho mais próximo.
lowess no R (ver Cleveland, 1979).Código da figura anterior.
A camada geom_smooth() tem como default o método lowess
ggplot(df_covid,aes(x=`DATA OBITO`,y=Mortes))+ geom_point(col="darkblue") + geom_smooth(col="red")
Podemos usar o argumento span para controlar a suavização.
Números menores produzem linhas mais onduladas, números maiores produzem linhas mais suaves.
ggplot(df_covid,aes(x=`DATA OBITO`,y=Mortes))+ geom_point( col="darkblue") + geom_smooth(col="red")+ geom_smooth(span=0.08, col="yellow", size=1.1)
ggplot(AIDS,aes(x=time(AIDS),y=as.matrix(AIDS)))+ geom_point(col="darkblue") + geom_smooth(col="red")
Técnicas de suavização podem ser aplicadas também no alisamento de uma série temporal em função de outra.
Aqui, vamos suavizar o diagrama de dispersão de duas séries temporais medidas simultaneamente.
Vamos considerar os dados de possíveis efeitos da temperatura (temp) sobre a mortalidade cardiovascular (cmort) em Los Angeles.
A figura seguinte mostra o diagrama de dispersão de mortalidade, \(M_t\), vs. temperatura, \(T_t\), e a linha suavizada via lowess.
Figura 16: Alisamento da mortalidade como função da temperatura via lowess.
Notem uma relação não-linear entre mortalidade e temperatura.
A mortalidade aumenta a temperaturas extremas, mas de uma forma assimétrica; é mais alta em temperaturas mais frias do que em temperaturas mais quentes.
A taxa de mortalidade mínima ocorre a aproximadamente 83\(^{\circ}\) F.
O código usado para reproduzir esta figura no R foi:
df_tm <- data.frame(tempr,cmort) ggplot(data = df_tm, aes(tempr, cmort)) + geom_point(color='blue') + labs(x = "\n Temperatura (em Fº) ", y = "Mortalidade \n") + geom_smooth(col="red")
Decomposição de séries é um procedimento matemático útil para dividir (decompor) a série em três componentes:
Tendência;
Sazonalidade;
Ruído (resíduo) aleatório.
Há dois tipos clássicos de decomposição; a aditiva,
série temporal = tendência + sazonalidade + ruído aleatório
e a multiplicativa,
série temporal = tendência * sazonalidade * ruído aleatório.
Para obter uma boa decomposição da série, o primeiro passo é escolher entre o modelo aditivo ou multiplicativo.
Duas regras básicas:
Se a variação sazonal parece constante ou não varia muito, aplicamos o modelo aditivo;
Se a magnitude da sazonalidade aumenta/diminui com o tempo, devemos aplicar o modelo multiplicativo.
Para ilustrar a técnica de decomposição, consideremos as duas séries temporais a seguir.
Figura 17: Produção de cerveja na Austrália.
A variação sazonal não varia com o tempo, i.e., é constante.
Então, vamos usar o modelo aditivo para esta série.
A variação sazonal aumenta ao longo do tempo.
Logo, devemos usar o modelo multiplicativo para esta série.
R útil para visualizar a variação sazonal de uma série é seasonplot do pacote forecast.seasonplot(st.ausbeer, s=4, col=rainbow(4)) seasonplot(a10, s=12, col=rainbow(12))
Este gráfico sazonal permite visualizar o padrão sazonal da série a cada ciclo, e.g., a cada ano, e é útil na identificação de anos em que o padrão muda.
Observem, por exemplo, o aumento na produção de cervejas na Austrália no 2º trimestre de cada ano.
Para as vendas de medicamentos, notem o aumento, principalmente no mês de janeiro.
Para ilustrar a técnica de decomposição, vamos usar algumas funções do pacote TSstudio e reproduzí-la passo-a-passo para entender como ela funciona.
Faremos o exemplo para demonstrar o passo-a-passo com a série que apresenta sazonalidade aditiva, st.ausbeer.
[Passo 1]
Inicialmente, vamos detectar a tendência das séries exemplificadas, via suavização média móvel.
Para decompor uma série, é crucial conhecer o período de sazonalidade: semanal, mensal, anual, etc.
A produção de cerveja australiana tem sazonalidade anual, ou seja, os ciclos parecem se repetir a cada ano.
Como é registrada trimestralmente, há 4 observações por ano, então, usamos janelas (vizinhos imediatos) de média móvel de 4.
[Passo 1]
A regra geral é usar média móvel de dois lados com ordem de frequência/2.
Por exemplo, no caso do dataset st.ausbeer, como a frequência é trimestral, vamos usar média móvel de dois lados, onde cada lado terá 2 observações consecutivas.
Isso significa que, para calcular o valor suavizado da observação \(t\), calcularemos a média da observação \(t\) junto com as 2 observações anteriores e seguintes (ou seja, com média de 5 pontos dos dados).
ts_ma {TSstudio} do R:ap_smooth <- ts_ma(st.ausbeer, n = 2, plot = FALSE,
separate = FALSE)
ap_smooth$plot %>%
layout(legend = list(x = 0.1, y = 0.9))
Na próxima etapa, estimamos a tendência
Converteremos a série e a tendência suavizada em um objeto data.frame com a função ts_to_prophet para mesclar as duas séries:
# modelo aditivo
ap_smooth <- ts_ma(st.ausbeer, n = 2, plot = FALSE, separate = FALSE)
df <- ts_to_prophet(st.ausbeer) %>%
select(date = ds, y) %>%
left_join(ts_to_prophet(ap_smooth$ma_2) %>%
select(date = ds, trend = y), by = "date")
df$detrend <- df$y - df$trend
df$year <- year(df$date)
df$month <- month(df$date)
ts_plot(df[, c("date", "y", "trend")], title = "Beer production and trend component",
type = "multiple")
[Passo 2]
# modelo aditivo
df$detrend <- df$y - df$trend
df$year <- year(df$date)
df$month <- month(df$date)
p <- plot_ly()
for(i in unique(df$year)){
temp <- NULL
temp <- df %>% filter(year == i)
p <- p %>% add_lines(x = temp$month,
y = temp$detrend,
name = i)
}
seasonal_comp <- df %>%
group_by(month) %>%
summarise(month_avg = mean(detrend, na.rm = TRUE),
.groups = "drop")
p %>% add_lines(x = seasonal_comp$month,
y = seasonal_comp$month_avg,
line = list(color = "black", dash = "dash", width = 4),
name = "Seasonal Component")
[Passo 3]
Vamos, então, obter a sazonalidade média da série.
Observem que em cada ano há um ciclo diferente na série das sazonalidades.
Cria-se uma matriz em que as linhas representam os anos e as colunas contêm os elementos do mesmo período (mesmo mês, mesmo trimestre, etc) e calcula a média de cada coluna.
Este passo 3 corresponde a traçar a curva sazonal média.
# modelo aditivo
seasonal_comp <- df %>%
group_by(month) %>%
summarise(month_avg = mean(detrend, na.rm = TRUE),
.groups = "drop")
df <- df %>% left_join(seasonal_comp, by = "month")
df$irregular <- df$y - df$trend - df$month_avg
ts_plot(df[, c("date", "month_avg")],
title = "Seasonal Component",
type = "multiple")
[Passo 4] Computando o ruído aleatório.
As etapas anteriores já extraíram muitas informações da série temporal original. Resta apenas o que chamamos de ruído aleatório.
Pela formulação do modelo, sabemos que na decomposição aditiva,
ruído aleatório = série temporal - tendência - sazonalidade
e na multiplicativa,
ruído aleatório = série temporal / (tendência * sazonalidade).
Ou seja, o ruído aleatório é o que resta de informação após eliminação da sazonalidade e da tendência da série temporal.
df <- df %>% left_join(seasonal_comp, by = "month") df$irregular <- df$y - df$trend - df$month_avg head(df)
ts_plot(df[, c("date", "y" ,"trend", "detrend", "month_avg", "irregular")],
title = "Beer Production and its Components",
type = "multiple")
Para tornar nossa vida mais fácil, alguns pacotes do R fornecem a decomposição com uma única linha de comando.
Por exemplo decompose, stl, STL.
As etapas apresentadas aqui reproduzem a função ts_decompose do pacote TSstudio.
Basta indicar o período correspondente ao ciclo sazonal da série e o modelo a ser aplicado.
ts_decompose(st.ausbeer, type = "additive") ts_decompose(a10, type = "multiplcative") ts_covid <- stats::ts(df_covid$Mortes, start=2020-03-28, frequency=365) ts_decompose(ts_covid, type = "both")
É uma medida de quão “fácil” a série é de prever.
De acordo com Hyndman e Athanasopoulos (2021), uma série com alta sazonalidade e tendência (bem comportada) terá uma entropia perto de 0, enquanto uma série com muitos ruídos, logo difícil de prever, terá uma entropia perto de 1.
Para obter essa medida basta usar a função feat_spectral do pacote feasts.
feat_spectral(st.ausbeer)
## spectral_entropy ## 0.3754978
feat_spectral(a10)
## spectral_entropy ## 0.3499013
feat_spectral(ts_covid)
## spectral_entropy ## 0.4194675
feat_spectral(cmort)
## spectral_entropy ## 0.7684916
tidyverse, ggplot2, gridExtra, lubridate, dplyr, reshape2, zoo, readxl, astsa, fpp, forecast, magrittr, rlang, tidyverse, lubridate, feasts, zoo, fpp, TSstudio, plotly.Acesse: https://github.com/rladiessalvador
Siga: @rladiessalvador